# Meta-analysis for healthy control (HC) vs ENL

# Datasets: GSE129033 and GSE74481

library("MetaDE")
library("readxl")
library("org.Hs.eg.db")
library("AnnotationDbi")
library("WriteXLS")

# Reading Expression Matrix from GSE129033

GSE129033 <- read.csv(".../GSE129033/GSE129033_healthy_vs_ENL.csv", header = TRUE)
head(GSE129033, 2)

# Formatting

table(duplicated(GSE129033$ENTREZID))
row.names(GSE129033) <- as.character(GSE129033$ENTREZID)
GSE129033 <- GSE129033[,c(13:24)]
head(GSE129033, 2)
label_GSE129033 <- ifelse(grepl("healthy", colnames(GSE129033)) == TRUE, '1', '0')
label_GSE129033

# Reading Expression Matrix from GSE74481

GSE74481 <- read.csv(".../GSE74481/GSE74481_healthy_vs_ENL.csv", header = TRUE)
head(GSE74481, 2)

# Formatting 

table(duplicated(GSE74481$ENTREZID))
row.names(GSE74481) <- as.character(GSE74481$ENTREZID)
GSE74481 <- GSE74481[, c(2:18)]
head(GSE74481, 2)
label_74481 <- ifelse(grepl("LL", colnames(GSE74481)) == TRUE, '1', '0') 
label_74481

#Constructing list

i <- Reduce(intersect, x = list(row.names(GSE129033), rownames(GSE74481)))

# exclude genes from ribosomal proteins
ribosome <- read.table("../ribosomal_proteins.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)
head(ribosome, 3)
i <- setdiff(i, ribosome$gene_id)

healthy_vs_ENL <- list(data.matrix(GSE129033[i, ]), data.matrix(GSE74481[i, ]))
rm(i)

names(healthy_vs_ENL) <- c("GSE129033", "GSE74481")
K <- length(healthy_vs_ENL)
label <- list(GSE129033 = label_129033, GSE74481 = label_74481)

clin.data <- lapply(label, function(x) {data.frame(x)} )
for (k in 1:length(clin.data)){
        colnames(clin.data[[k]]) <- "label"
}

# meta analysis using REM
meta.res <- MetaDE(data = healthy_vs_ENL, clin.data = clin.data, data.type = "continuous", resp.type = "twoclass", response = 'label', ind.method = rep('limma', 4), meta.method = "REM", select.group = c('0', '1'), ref.level = c('0'), paired = rep(FALSE, length(healthy_vs_ENL)), REM.type = "HO", tail = 'abs')

# saving main holder object 

saveRDS(meta.res, file = "obj/healthy_vs_ENL_meta.rds", compress = "bzip2")

# Constructing a matrix to store results

healthy_vs_ENL.rem2 <- matrix(data = c(meta.res$meta.analysis$mu.hat, meta.res$meta.analysis$mu.var, meta.res$meta.analysis$tau2, meta.res$meta.analysis$FDR), nrow = length(meta.res$meta.analysis$mu.hat), ncol = 4, dimnames = list(row.names(meta.res$meta.analysis$FDR), c("muhat", "muvar", "tau2", "FDR")), byrow = FALSE)
healthy_vs_ENL.rem2 <- cbind(meta.res$ind.ES[row.names(healthy_vs_ENL.rem2), ], healthy_vs_ENL.rem2)
healthy_vs_ENL.rem2 <- as.data.frame(healthy_vs_ENL.rem2)
healthy_vs_ENL.rem2$Symbol <- select(org.Hs.eg.db, keys = row.names(healthy_vs_ENL.rem2),columns = c("SYMBOL"), keytype = "ENTREZID")$SYMBOL 
colnames(healthy_vs_ENL.rem2)[1:4] <- names(healthy_vs_ENL)
healthy_vs_ENL.rem2 <- merge(healthy_vs_ENL.rem2, meta.res$ind.Var, by = "row.names")
colnames(healthy_vs_ENL.rem2)[11:14] <- paste("var", names(healthy_vs_ENL), sep = "_")

WriteXLS(healthy_vs_ENL.rem2, ExcelFileName = "results/healthy_vs_ENL_ES.xls")

rm(list = setdiff(ls(), c("meta.res", 'draw.DEnumber_all')))

# Generating P-values using limma

p_GSE129033 <- read.csv("GSE129033_healthy_vs_ENL.csv", header = TRUE)
p_GSE129033 <- p_GSE129033[,c("ENTREZID", "P_value", "adjusted_P_value")]
p_GSE129033 <- p_GSE129033[!is.na(p_GSE129033$ENTREZID),]
row.names(p_GSE129033) <- p_GSE129033$ENTREZID

p_GSE74481 <- as.data.frame(read_xls("../../GSE74481/results/results.healthy_vs_ENL.xls"))
p_GSE74481 <- p_GSE74481[,c("ENTREZID", "P_value", "adjusted_P_value")]
p_GSE74481 <- p_GSE74481[!is.na(p_GSE74481$ENTREZID),]
row.names(p_GSE74481) <- p_GSE74481$ENTREZID

## Genes common to both the studies
i <- Reduce(intersect, list(row.names(p_GSE129033), row.names(p_GSE74481)))

## exclude genes from ribosomal proteins
ribosome <- read.table("../ribosomal_proteins.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)
i <- setdiff(i, ribosome$gene_id)
tmp <- cbind(p_GSE129033[i,]$P_value, p_GSE74481[i,]$P_value)
row.names(tmp) <- i
colnames(tmp) <- c("GSE129033", "GSE74481")
rm(i)

# Storing P_values from both studies
healthy_vs_ENL_pvalues <- list(p = tmp)
rm(tmp)

# Meta-analysis using other methods

healthy_vs_ENL<- MetaDE.pvalue(healthy_vs_ENL_pvalues, meta.method = c("roP","maxP", "SR", "fisher"), rth = 4, parametric = FALSE)

#Merging results

res <- cbind(healthy_vs_ENL$ind.p, healthy_vs_ENL$meta.analysis$pval)
tmp <- cbind(meta.res$meta.analysis$pval, meta.res$meta.analysis$FDR)
res <- cbind(res, tmp[row.names(res),][,1])
colnames(res)[8] <- "REM" 
rm(tmp)

i <- Reduce(intersect, list(row.names(p_GSE129033), row.names(p_GSE74481)))
tmp <- cbind(p_GSE129033[i,]$adjusted_P_val, p_GSE74481[i,]$adjusted_P_val)
row.names(tmp) <- i
colnames(tmp) <- c("GSE129033","GSE74481")
res <- cbind(tmp[i, ], healthy_vs_ENL$meta.analysis$FDR[i, ])
res <- cbind(res[i, ], meta.res$meta.analysis$FDR[i,])
colnames(res)[8] <- "REM" 
rm(tmp)
res <- as.data.frame(res)
res$Symbol <- select(org.Hs.eg.db, keys = row.names(res), columns = "SYMBOL", "ENTREZID")$SYMBOL 

# Saving results from other methods
WriteXLS(res, ExcelFileName = ".../healthy_vs_ENL_allMethods.xls", row.names = TRUE)


# Intersection between all methods 
intersection <- data.frame(ENTREZID = Reduce(intersect, list(sdef, REM, roP, maxP, SR, fisher)))
intersection$SYMBOL <- select(org.Hs.eg.db, keys = as.character(intersection$ENTREZID), columns = "SYMBOL", "ENTREZID")$SYMBOL
WriteXLS(intersection, ExcelFileName = ".../Intersection.xls")

# saving the intersection between all methods except REM
intersection <- data.frame(ENTREZID = Reduce(intersect, list(sdef, SR, fisher, maxP)))
intersection$SYMBOL <- select(org.Hs.eg.db, keys = as.character(intersection$ENTREZID), columns = "SYMBOL", "ENTREZID")$SYMBOL
WriteXLS(intersection, ExcelFileName = ".../Intersection_except_REM.xls")